library(gtsummary)
library(dplyr)
library(tidyverse)
library(survival)
library(survminer)
library(lubridate)
library(readxl)
library(janitor)
library(ggplot2)
library(tidycmprsk)
library(ggsurvfit)
library(tidymodels)
library(sjPlot)
library(patchwork)
library(bstfun)
library(ComplexUpset)
library(writexl)
library(rms)
library(timeROC)
library(rms)
library(pec)
library(prodlim)
library(riskRegression)
library(plotly)
library(ggtext)
library(WeightIt)
library(cobalt)
library(survey)
db_Kansas = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data Template 2.11.2025 - due 5-1-25 (de-identified).xlsx", sheet = "bsAb_DATA")%>%
filter(`bsAb name` %in% c("1","3"))%>%
clean_names()%>%
filter(!is.na(age_at_bs_ab_initiation)) %>%
mutate(
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
baseline_anc = as.double(baseline_anc),
baseline_alc = as.double(baseline_alc),
alc_day_8 = as.double(alc_day_8),
alc_day_15 = as.double(alc_day_15),
alc_day_22 = as.double(alc_day_22),
alc_c2d1 = as.double(alc_c2d1),
anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
)
# Write to Excel
#write_xlsx(db_Kansas, "Kansas_elra_data.xlsx")
db_CC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/CC bsAb Consortium Submission CCF 6-11-25.xlsx", sheet = "bsAb_DATA")%>%
filter(`bsAb name` %in% c("1","3"))%>%
clean_names()%>%
mutate(
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
baseline_anc = as.double(baseline_anc),
baseline_alc = as.double(baseline_alc),
alc_day_8 = as.double(alc_day_8),
alc_day_15 = as.double(alc_day_15),
alc_day_22 = as.double(alc_day_22),
alc_c2d1 = as.double(alc_c2d1),
anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
)
# Write to Excel
#write_xlsx(db_CC, "CC_elra_data.xlsx")
db_HUMC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/BsAb Consortium Data HUMC.xlsx", sheet = "bsAb_DATA")%>%
filter(`bsAb name` %in% c("1","3"))%>%
clean_names()%>%
mutate(
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
baseline_anc = as.double(baseline_anc),
baseline_alc = as.double(baseline_alc),
alc_day_8 = as.double(alc_day_8),
alc_day_15 = as.double(alc_day_15),
alc_day_22 = as.double(alc_day_22),
alc_c2d1 = as.double(alc_c2d1),
anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
)
# Write to Excel
#write_xlsx(db_HUMC, "HUMC_elra_data.xlsx")
db_MCC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data MCC-deidentified.xlsx", sheet = "bsAb_DATA")%>%
filter(`bsAb name` %in% c("1","3"))%>%
clean_names()%>%
mutate(
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
baseline_anc = as.double(baseline_anc),
baseline_alc = as.double(baseline_alc),
cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
)
db_Huntsman = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/Utah_Huntsman_bsAb Consortium Data FINAL_5.4.25.xlsx", sheet = "bsAb_DATA")%>%
filter(`bsAb name` %in% c("1","3"))%>%
clean_names()%>%
mutate(
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
baseline_anc = as.double(baseline_anc)/1000,
baseline_alc = as.double(baseline_alc)/1000,
anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
)
# Write to Excel
#write_xlsx(db_Huntsman, "Huntsman_elra_data.xlsx")
db_FHCC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/FHCC Bispecifics (5.5.2025).xlsx", sheet = "bsAb_DATA")%>%
filter(`bsAb name` %in% c("1","3"))%>%
clean_names()%>%
mutate(
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
date_of_infection_80 = ymd(date_of_infection_80),
baseline_anc = as.double(baseline_anc),
cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
)
db_UTSW = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data Template 2.11.2025 1.xlsx", sheet = "bsAb_DATA")%>%
filter(`bsAb name` %in% c("1","3")) %>%
clean_names()%>%
mutate(
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
baseline_anc = as.double(baseline_anc),
baseline_alc = as.double(baseline_alc),
anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles)
)
# Write to Excel
#write_xlsx(db_UTSW, "UTSW_elra_data.xlsx")
db_MDACC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/MDACC BITE ASH 2025.xlsx", sheet = "bsAb_DATA")%>%
filter(`bsAb name` %in% c("1","3"))%>%
clean_names()%>%
mutate(
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
baseline_anc = as.double(baseline_anc),
baseline_alc = as.double(baseline_alc),
anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1),
total_number_paraskelatal_plasmacytomas = as.double(total_number_paraskelatal_plasmacytomas)
)
# Write to Excel
#write_xlsx(db_MDACC, "MDACC_elra_data.xlsx")
db_Roswell = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data Template 5.8.2025_NN.xlsx", sheet = "bsAb_DATA")%>%
filter(`bsAb name` %in% c("1","3"))%>%
clean_names()%>%
mutate(
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
baseline_anc = as.double(baseline_anc),
baseline_alc = as.double(baseline_alc),
anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1),
total_number_paraskelatal_plasmacytomas = as.double(total_number_paraskelatal_plasmacytomas)
)
# Write to Excel
#write_xlsx(db_Roswell, "Roswell_elra_data.xlsx")
# Combine all into a single dataframe
combined_db <- bind_rows(
list(
HUMC = db_HUMC,
MCC = db_MCC,
Huntsman = db_Huntsman,
FHCC = db_FHCC,
UTSW = db_UTSW,
MDACC = db_MDACC,
Roswell = db_Roswell,
CC = db_CC,
Kansas = db_Kansas
),
.id = "site"
) %>%
mutate(
prior_bcma_directed_therapy_yes_no = ifelse(!is.na(most_recent_bcma_agent_name),1,0),
bs_ab_name = case_when(
bs_ab_name == 1 ~ "tec",
bs_ab_name == 3 ~ "elra",
.default = NULL
),
max_crs_grade_0_5 = ifelse(is.na(max_crs_grade_0_5), 0,max_crs_grade_0_5), ## Set NAs to 0
max_icans_grade = ifelse(is.na(max_icans_grade), 0, max_icans_grade ) ## Set NAs to 0
)
# Convert all datetime columns to Date
combined_db_clean <- combined_db %>%
mutate(across(where(lubridate::is.POSIXt), as.Date))
# Write to Excel
write_xlsx(combined_db_clean, "elra_tec_data.xlsx")
# Write to Excel
# write_xlsx(combined_db_clean %>% filter(prior_bcma_directed_therapy_yes_no != prior_bcma_directed_therapy), "elra_data_discordant.xlsx")
## Several institutions do not have any cases of elranatamab
# db_UABMC = read_excel(path = "Data/bsAb Consortium Data 5.8.xlsx", sheet = "bsAb_DATA")%>%
# filter(`bsAb name` == "3")%>%
# clean_names()
# db_Stanford = read_excel(path = "Data/De-identified Stanford bsAb Consortium Data Template 4.30.2025.xlsx", sheet = "bsAb_DATA")%>%
# filter(`bsAb name` == "3")%>%
# clean_names()
# db_LCI = read_excel(path = "Data/bsAb Consortium Data LCI.xlsx", sheet = "bsAb_DATA")%>%
# filter(`bsAb name` == "3")%>%
# clean_names()
# db_COH = read_excel(path = "Data/bsAb Consortium Data Template 2.11.2025_COHv2.xlsx", sheet = "bsAb_DATA")%>%
# filter(`bsAb name` == "3")%>%
# clean_names()
# db_MUSC = read_excel(path = "Data/USMM bsAb Consortium Data Template 2.11.2025.xlsx", sheet = "bsAb_DATA")%>%
# filter(`bsAb name` == "3")%>%
# clean_names()
# db_Mayo = read_excel(path = "Data/Bispecific_consortium_20250521 NoPHI.xlsx", sheet = "bsAb_DATA") %>%
# filter(`bsAb name` == 3)%>%
# clean_names()
combined_db <- combined_db %>%
mutate(
# baseline_anc = as.numeric(baseline_anc),
# baseline_hgb,
# baseline_pl_ts,
# baseline_crp_prior_to_bs_ab_initiation,
# baseline_ferritin_prior_to_bs_ab_initiation,
HT_Score =
(baseline_anc <= 1.2) +
(baseline_hgb <= 9.0) +
(baseline_pl_ts >= 76 & baseline_pl_ts <= 175) +
(baseline_crp_prior_to_bs_ab_initiation >= 3.0) +
(baseline_ferritin_prior_to_bs_ab_initiation >= 650 & baseline_ferritin_prior_to_bs_ab_initiation <= 2000) +
(baseline_pl_ts <= 75) + ## Add an extra point for Plt <=75
(baseline_ferritin_prior_to_bs_ab_initiation > 2000), ## Add an extra point for ferritin >2000
# Classify the risk based on the score
HT_Risk = factor(
case_when(
HT_Score >= 2 ~ "HThigh",
baseline_pl_ts <=75 ~ "HThigh", ## Equivalent to at least 2 points
baseline_ferritin_prior_to_bs_ab_initiation > 2000 ~ "HThigh", ## Equivalent to at least 2 points
(baseline_anc <= 1.2) + (baseline_hgb <= 9.0) + (baseline_pl_ts >= 76 & baseline_pl_ts <= 175) >=2 ~ "HThigh",
HT_Score <2 ~ "HTlow",
.default = "NA"
),
levels = c("HTlow", "HThigh")
)
)
write_xlsx(combined_db %>% select(site, patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate, baseline_anc, baseline_hgb, baseline_pl_ts, baseline_crp_prior_to_bs_ab_initiation, baseline_ferritin_prior_to_bs_ab_initiation, HT_Score, HT_Risk), "HT_score.xlsx")
combined_db <- combined_db %>%
mutate(
severe_infection = ifelse(infection_severity_73 %in% 1 | infection_severity_76 %in% 1 | infection_severity_79 %in% 1| infection_severity_82 %in% 1,1,0),
first_severe_infection_date = case_when(
infection_severity_73 %in% 1 ~ date_of_infection_74,
infection_severity_76 %in% 1 ~ date_of_infection_77,
infection_severity_79 %in% 1 ~ date_of_infection_80,
infection_severity_82 %in% 1 ~ date_of_infection_83,
.default = NULL
),
first_severe_infection_type = case_when(
infection_severity_73 %in% 1 ~ number_infection_1_bacterial_viral_fungal,
infection_severity_76 %in% 1 ~ number_infection_2_bacterial_viral_fungal,
infection_severity_79 %in% 1 ~ number_infection_3_bacterial_viral_fungal,
infection_severity_82 %in% 1 ~ number_infection_4_bacterial_viral_fungal,
.default = NULL
)
)
theme_gtsummary_compact()
## Setting theme "Compact"
variables <- c("elra", "Age", "Female sex", "ECOG 2+", "EMD", "Hgb", "LDH", "Prior BCMA")
data1 <- combined_db %>%
drop_na(age_at_bs_ab_initiation, ecog_ps_at_bs_ab_initiation_0_5, just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk, baseline_hgb, baseline_ldh_prior_to_bs_ab_initiation) %>%
transmute(
Age = age_at_bs_ab_initiation,
"Female sex" = sex_0_male_1_female,
elra = ifelse(bs_ab_name == "elra", 1, 0),
"ECOG 2+" = ifelse(ecog_ps_at_bs_ab_initiation_0_5 >=2,1,0),
EMD = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,1,0),
Hgb = as.double(baseline_hgb),
LDH = as.double(baseline_ldh_prior_to_bs_ab_initiation),
"Prior BCMA" = prior_bcma_directed_therapy_yes_no
)
uv_tab <- tbl_uvregression(
data1[c(variables)],
method = glm,
y = "elra",
method.args = list(family = binomial),
exponentiate = TRUE
) %>%
sort_p()
mv_tab<-glm(elra ~ Age + `Female sex` + `ECOG 2+` + EMD + Hgb + LDH + `Prior BCMA`, data=data1, family=binomial) %>%
tbl_regression(exponentiate=TRUE)
tbl_merge(list(uv_tab, mv_tab), tab_spanner = c("**Univariate**", "**Multivariable**"))
| Characteristic |
Univariate
|
Multivariable
|
|||||
|---|---|---|---|---|---|---|---|
| N | OR1 | 95% CI1 | p-value | OR1 | 95% CI1 | p-value | |
| Hgb | 371 | 1.12 | 1.01, 1.25 | 0.029 | 1.14 | 1.02, 1.28 | 0.026 |
| EMD | 371 | 1.69 | 0.95, 2.96 | 0.069 | 1.71 | 0.96, 3.04 | 0.068 |
| Age | 371 | 1.02 | 1.00, 1.04 | 0.087 | 1.01 | 0.99, 1.04 | 0.3 |
| Female sex | 371 | 1.43 | 0.93, 2.21 | 0.11 | 1.43 | 0.92, 2.24 | 0.11 |
| Prior BCMA | 371 | 0.76 | 0.49, 1.17 | 0.2 | 0.78 | 0.50, 1.23 | 0.3 |
| ECOG 2+ | 371 | 1.19 | 0.75, 1.88 | 0.4 | 1.27 | 0.78, 2.07 | 0.3 |
| LDH | 371 | 1.00 | 1.00, 1.00 | 0.8 | 1.00 | 1.00, 1.00 | >0.9 |
| 1 OR = Odds Ratio, CI = Confidence Interval | |||||||
IPTW.data <- combined_db %>%
drop_na(age_at_bs_ab_initiation, ecog_ps_at_bs_ab_initiation_0_5, just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk, baseline_hgb, baseline_ldh_prior_to_bs_ab_initiation) %>%
mutate(
Age = age_at_bs_ab_initiation,
"Female sex" = sex_0_male_1_female,
bsAb = bs_ab_name,
"Received full dose" = ifelse(!is.na(date_of_first_bs_ab_full_dose),1,0),
"ECOG" = case_when(
ecog_ps_at_bs_ab_initiation_0_5 %in% c("0","1") ~"0-1",
ecog_ps_at_bs_ab_initiation_0_5 ==2 ~"2",
ecog_ps_at_bs_ab_initiation_0_5 >=3 ~"3+",
.default = NULL),
EMD = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,1,0),
Hgb = as.double(baseline_hgb),
LDH = as.double(baseline_ldh_prior_to_bs_ab_initiation),
"Prior BCMA" = prior_bcma_directed_therapy_yes_no
)
## Use the WeightIt function to generate propensity scores and weights
W.out <- weightit(bsAb ~ Age + `ECOG` + EMD + Hgb + LDH + `Prior BCMA`,
data=IPTW.data,
method = "glm", ## Methods: glm, gbm, cbps, npcbps, ebal, optweight, super, bart, energy
estimand = "ATE",
stabilize = TRUE)
bal.tab(W.out, stats = c("m", "v"), thresholds = c(m=0.25), binary = "std")
## Balance Measures
## Type Diff.Adj M.Threshold V.Ratio.Adj
## prop.score Distance 0.0050 Balanced, <0.25 0.9733
## Age Contin. 0.0090 Balanced, <0.25 0.9595
## ECOG_0-1 Binary -0.0121 Balanced, <0.25 .
## ECOG_2 Binary 0.0152 Balanced, <0.25 .
## ECOG_3+ Binary -0.0027 Balanced, <0.25 .
## EMD Binary -0.0034 Balanced, <0.25 .
## Hgb Contin. -0.0186 Balanced, <0.25 1.0974
## LDH Contin. -0.0276 Balanced, <0.25 1.0961
## Prior BCMA Binary -0.0119 Balanced, <0.25 .
##
## Balance tally for mean differences
## count
## Balanced, <0.25 9
## Not Balanced, >0.25 0
##
## Variable with the greatest mean difference
## Variable Diff.Adj M.Threshold
## LDH -0.0276 Balanced, <0.25
##
## Effective sample sizes
## elra tec
## Unadjusted 123. 248.
## Adjusted 113.38 242.39
summary(W.out)
## Summary of weights
##
## - Weight ranges:
##
## Min Max
## elra 0.4745 |---------------------------| 2.3134
## tec 0.7811 |-------------| 1.7418
##
## - Units with the 5 most extreme weights by group:
##
## 123 102 81 91 9
## elra 1.5458 1.5896 1.653 1.7607 2.3134
## 169 71 274 258 271
## tec 1.4257 1.4453 1.5215 1.6783 1.7418
##
## - Weight statistics:
##
## Coef of Var MAD Entropy # Zeros
## elra 0.292 0.222 0.040 0
## tec 0.152 0.115 0.011 0
##
## - Effective Sample Sizes:
##
## elra tec
## Unweighted 123. 248.
## Weighted 113.38 242.39
love_plot <- love.plot(W.out,
binary = "std",
line = TRUE,
abs = TRUE,
thresholds = c(m=0.2),
sample.names = c("Unweighted", "IPTW"),
title = "A)",
var.order = "unadjusted",
var.names = c(
prop.score = "Distance"
))+ theme(plot.title = element_text(hjust = 0))
love_plot
weight_df = data.frame(W.out$weights, W.out$treat)
IPTW.data$weights <- W.out$weights
IPTW.data$ps <- W.out$ps
write.csv(data1, "IPTW.dataset.csv")
Weights_plot <- ggplot(
IPTW.data, aes(x=ps, y = weights, color = bsAb)) +
geom_point(shape = 16) +
theme_classic() +
scale_y_continuous(breaks = seq(0,4 , by = 0.5)) +
labs(x = "Propensity score", y = "Weights") +
ggtitle("B)")
Weights_plot
theme_gtsummary_compact()
combined_db %>%
transmute(
"Site" = site,
bs_ab_name,
) %>%
tbl_summary(
by = bs_ab_name,
missing = "no",
sort = list(all_categorical() ~ "frequency"),
statistic = list(
all_continuous() ~ "{median} ({p25}, {p75})",
all_categorical() ~ "{n} ({p}%)"
)
) %>%
bold_labels()
| Characteristic | elra N = 1301 |
tec N = 2771 |
|---|---|---|
| Site | ||
| Â Â Â Â MCC | 53 (41%) | 72 (26%) |
| Â Â Â Â MDACC | 12 (9.2%) | 57 (21%) |
| Â Â Â Â CC | 20 (15%) | 47 (17%) |
| Â Â Â Â UTSW | 10 (7.7%) | 40 (14%) |
| Â Â Â Â FHCC | 7 (5.4%) | 24 (8.7%) |
| Â Â Â Â Huntsman | 9 (6.9%) | 22 (7.9%) |
| Â Â Â Â Kansas | 4 (3.1%) | 13 (4.7%) |
| Â Â Â Â HUMC | 10 (7.7%) | 1 (0.4%) |
| Â Â Â Â Roswell | 5 (3.8%) | 1 (0.4%) |
| 1 n (%) | ||
theme_gtsummary_compact()
combined_db %>%
transmute(
#"Site" = site,
"Age (years)" = age_at_bs_ab_initiation,
"Female sex" = sex_0_male_1_female == 1,
"Race" = factor(case_when(
is.na(race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5) ~ NA,
race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 0 ~ "White",
race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 1 ~ "Black",
race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 2 ~ "AAPI",
race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 3 ~ "American Indian",
race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 4 ~ "Other",
race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 5 ~ NA
), levels = c("White","Black","AAPI","American Indian","Other",NA)),
#ecog_ps_at_bs_ab_initiation_0_5,
"ECOG" = case_when(
ecog_ps_at_bs_ab_initiation_0_5 == 0 ~ "0",
ecog_ps_at_bs_ab_initiation_0_5 == 1 ~ "1",
ecog_ps_at_bs_ab_initiation_0_5 == 2 ~ "2",
ecog_ps_at_bs_ab_initiation_0_5 >= 3 ~ "3+"
),
"LDH (U/L)" = as.double(baseline_ldh_prior_to_bs_ab_initiation),
"Hgb (g/dL)" = baseline_hgb,
"Heavy chain" = factor(case_when(
myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgG" ~ "IgG",
myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgM" ~ "IgM",
myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgA" ~ "IgA",
myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgD" ~ "IgD",
.default = "None"
), levels = c("IgG", "IgA","IgM","IgD","None")),
"True EMD" = just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,
"Oligo or non-secretory" = oligo_or_non_secretory_at_bs_ab_initiation_0_no_1_yes,
"≥1 HRCA" = ifelse(deletion_17p_prior_to_bs_ab_0_no_1_yes==1 | t_4_14_prior_to_bs_ab_0_no_1_yes == 1 | t_14_16_prior_to_bs_ab_0_no_1_yes == 1 | amp_1_q_only_0_no_1_yes == 1 | gain_amp1q21_prior_to_bs_ab_0_no_1_yes == 1,1,0),
"del(17p)" = deletion_17p_prior_to_bs_ab_0_no_1_yes,
"t(4;14)" = t_4_14_prior_to_bs_ab_0_no_1_yes,
"t(14;16)" = t_14_16_prior_to_bs_ab_0_no_1_yes,
"gain/amp(1q)" = gain_amp1q21_prior_to_bs_ab_0_no_1_yes,
"amp(1q)" = amp_1_q_only_0_no_1_yes,
"Prior LOTs" = number_of_prior_lines_of_therapy_before_bs_ab,
"Prior ASCT" = prior_auto_sct_no_0_yes_1,
"Prior CAR-T" = prior_cart_0_no_1_yes,
"Bortezomib refractory" = bortezomib_status == 1,
"Lenalidomide refractory" = lenalidomide_status == 1,
"Pomalidomide refractory" = pomalidomide_status == 1,
"Anti-CD38 refractory" = daratumumab_or_isatuximab_status == 1,
"Triple refractory" = triple_class_refractory_i_mi_d_pi_anti_cd38 == 1,
"Penta refractory" = penta_refractory_2_p_is_2_i_mi_ds_anti_cd38 == 1,
"Prior BCMA-directed therapy" = !is.na(most_recent_bcma_agent_name),
"Most recent BCMA" = factor(case_when(
most_recent_bcma_agent_name %in% c("Abecma", "ABECMA","Abecma on trial","Ide-cel") ~ "CAR-T",
most_recent_bcma_agent_name %in% c("Carvykti","carvykti", "Carvikty", "Clita-cel", "Cilta-cel") ~ "CAR-T",
most_recent_bcma_agent_name %in% c("CC98633 CAR-T","ddBCMA") ~ "CAR-T",
most_recent_bcma_agent_name %in% c("BLENREP", "blenrep", "Belantamab") ~ "ADC",
most_recent_bcma_agent_name %in% c("Teclistamab") ~ "TCE"
), levels = c("CAR-T","TCE", "ADC")),
"Trial eligible" = eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you == 1,
"Received full dose" = ifelse(!is.na(date_of_first_bs_ab_full_dose),1,0),
"Baseline ferritin (ng/mL)" = as.double(baseline_ferritin_prior_to_bs_ab_initiation),
"Baseline Hgb" = as.double(baseline_hgb),
"Baseline LDH" = as.double(baseline_ldh_prior_to_bs_ab_initiation),
"Duration of prior BCMA" = as.numeric(ymd(date_of_last_exposure_to_prior_bcma_agent) - ymd(date_of_first_exposure_to_prior_bcma_agent)),
bs_ab_name
) %>%
tbl_summary(
by = bs_ab_name,
missing = "no",
#sort = list(all_categorical() ~ "frequency"),
statistic = list(
all_continuous() ~ "{median} ({p25}, {p75})",
all_categorical() ~ "{n} ({p}%)",
"Age (years)" ~ "{median} ({min} to {max})"
)
) %>%
add_p() %>%
bold_labels() %>%
add_variable_grouping(
"High-risk cytogenetic abnormality at any time" = c("del(17p)","t(4;14)", "t(14;16)", "gain/amp(1q)", "amp(1q)"),
"Refractory status" = c("Bortezomib refractory", "Lenalidomide refractory", "Pomalidomide refractory", "Anti-CD38 refractory")
)
| Characteristic | elra N = 1301 |
tec N = 2771 |
p-value2 |
|---|---|---|---|
| Age (years) | 71 (39 to 95) | 69 (35 to 88) | 0.085 |
| Female sex | 74 (57%) | 128 (46%) | 0.044 |
| Race | 0.14 | ||
| Â Â Â Â White | 102 (82%) | 211 (76%) | |
| Â Â Â Â Black | 16 (13%) | 54 (20%) | |
| Â Â Â Â AAPI | 3 (2.4%) | 7 (2.5%) | |
| Â Â Â Â American Indian | 2 (1.6%) | 0 (0%) | |
| Â Â Â Â Other | 2 (1.6%) | 4 (1.4%) | |
| ECOG | 0.5 | ||
| Â Â Â Â 0 | 22 (17%) | 52 (20%) | |
| Â Â Â Â 1 | 60 (48%) | 130 (49%) | |
| Â Â Â Â 2 | 30 (24%) | 65 (25%) | |
| Â Â Â Â 3+ | 14 (11%) | 18 (6.8%) | |
| LDH (U/L) | 223 (177, 277) | 208 (166, 287) | 0.4 |
| Hgb (g/dL) | 10.05 (8.70, 11.70) | 9.60 (8.30, 11.30) | 0.020 |
| Heavy chain | <0.001 | ||
| Â Â Â Â IgG | 75 (58%) | 122 (44%) | |
| Â Â Â Â IgA | 28 (22%) | 47 (17%) | |
| Â Â Â Â IgM | 2 (1.5%) | 1 (0.4%) | |
| Â Â Â Â IgD | 2 (1.5%) | 1 (0.4%) | |
| Â Â Â Â None | 23 (18%) | 106 (38%) | |
| True EMD | 28 (22%) | 36 (13%) | 0.028 |
| Oligo or non-secretory | 0.3 | ||
| Â Â Â Â 0 | 116 (91%) | 259 (94%) | |
| Â Â Â Â 1 | 12 (9.4%) | 15 (5.4%) | |
| Â Â Â Â 2 | 0 (0%) | 2 (0.7%) | |
| ≥1 HRCA | 85 (69%) | 164 (63%) | 0.2 |
| High-risk cytogenetic abnormality at any time | |||
| Â Â Â Â del(17p) | 31 (25%) | 70 (27%) | 0.6 |
| Â Â Â Â t(4;14) | 17 (13%) | 34 (13%) | >0.9 |
| Â Â Â Â t(14;16) | 8 (6.3%) | 15 (5.8%) | 0.8 |
| Â Â Â Â gain/amp(1q) | 70 (56%) | 128 (49%) | 0.2 |
| Â Â Â Â amp(1q) | 9 (7.9%) | 18 (6.9%) | 0.7 |
| Prior LOTs | 6.00 (4.00, 8.00) | 6.00 (4.00, 8.00) | 0.8 |
| Prior ASCT | 66 (51%) | 163 (59%) | 0.13 |
| Prior CAR-T | 54 (42%) | 103 (37%) | 0.4 |
| Refractory status | |||
| Â Â Â Â Bortezomib refractory | 87 (67%) | 226 (82%) | <0.001 |
| Â Â Â Â Lenalidomide refractory | 110 (85%) | 238 (86%) | 0.7 |
| Â Â Â Â Pomalidomide refractory | 96 (74%) | 226 (82%) | 0.062 |
| Â Â Â Â Anti-CD38 refractory | 116 (89%) | 267 (97%) | 0.002 |
| Triple refractory | 118 (91%) | 253 (92%) | 0.8 |
| Penta refractory | 64 (49%) | 141 (51%) | 0.7 |
| Prior BCMA-directed therapy | 64 (49%) | 152 (55%) | 0.3 |
| Most recent BCMA | 0.012 | ||
| Â Â Â Â CAR-T | 54 (84%) | 57 (70%) | |
| Â Â Â Â TCE | 6 (9.4%) | 5 (6.2%) | |
| Â Â Â Â ADC | 4 (6.3%) | 19 (23%) | |
| Trial eligible | 28 (22%) | 63 (23%) | 0.8 |
| Received full dose | 121 (93%) | 268 (97%) | 0.093 |
| Baseline ferritin (ng/mL) | 293 (132, 832) | 393 (131, 1,254) | 0.3 |
| Baseline Hgb | 10.05 (8.70, 11.70) | 9.60 (8.30, 11.30) | 0.020 |
| Baseline LDH | 223 (177, 277) | 208 (166, 287) | 0.4 |
| Duration of prior BCMA | 0 (0, 22) | 0 (0, 126) | 0.034 |
| 1 Median (Min to Max); n (%); Median (Q1, Q3) | |||
| 2 Wilcoxon rank sum test; Pearson’s Chi-squared test; Fisher’s exact test | |||
theme_gtsummary_compact()
## Setting theme "Compact"
data1 = IPTW.data %>%
transmute(
"Age (years)" = age_at_bs_ab_initiation,
"Female sex" = sex_0_male_1_female == 1,
"Race" = factor(case_when(
is.na(race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5) ~ NA,
race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 0 ~ "White",
race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 1 ~ "Black",
race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 2 ~ "AAPI",
race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 3 ~ "American Indian",
race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 4 ~ "Other",
race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 5 ~ NA
), levels = c("White","Black","AAPI","American Indian","Other",NA)),
#ecog_ps_at_bs_ab_initiation_0_5,
"ECOG" = case_when(
ecog_ps_at_bs_ab_initiation_0_5 == 0 ~ "0",
ecog_ps_at_bs_ab_initiation_0_5 == 1 ~ "1",
ecog_ps_at_bs_ab_initiation_0_5 == 2 ~ "2",
ecog_ps_at_bs_ab_initiation_0_5 >= 3 ~ "3+"
),
"LDH (U/L)" = as.double(baseline_ldh_prior_to_bs_ab_initiation),
"Hgb (g/dL)" = baseline_hgb,
"Baseline ferritin (ng/mL)" = as.double(baseline_ferritin_prior_to_bs_ab_initiation),
"Heavy chain" = factor(case_when(
myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgG" ~ "IgG",
myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgM" ~ "IgM",
myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgA" ~ "IgA",
myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgD" ~ "IgD",
.default = "None"
), levels = c("IgG", "IgA","IgM","IgD","None")),
"True EMD" = just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,
"≥1 HRCA" = ifelse(deletion_17p_prior_to_bs_ab_0_no_1_yes==1 | t_4_14_prior_to_bs_ab_0_no_1_yes == 1 | t_14_16_prior_to_bs_ab_0_no_1_yes == 1 | amp_1_q_only_0_no_1_yes == 1 | gain_amp1q21_prior_to_bs_ab_0_no_1_yes == 1,1,0),
"Prior LOTs" = number_of_prior_lines_of_therapy_before_bs_ab,
"Prior ASCT" = prior_auto_sct_no_0_yes_1,
"Prior CAR-T" = prior_cart_0_no_1_yes,
"Triple refractory" = triple_class_refractory_i_mi_d_pi_anti_cd38 == 1,
"Penta refractory" = penta_refractory_2_p_is_2_i_mi_ds_anti_cd38 == 1,
"Prior BCMA-directed therapy" = !is.na(most_recent_bcma_agent_name),
"Trial eligible" = eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you == 1,
"Received full dose" = ifelse(!is.na(date_of_first_bs_ab_full_dose),1,0),
bs_ab_name,
weights
)
table_unweighted <- data1 %>% select(-weights) %>%
tbl_summary(
by = bs_ab_name,
missing = "no",
#sort = list(all_categorical() ~ "frequency"),
statistic = list(
all_continuous() ~ "{median} ({p25}, {p75})",
all_categorical() ~ "{n} ({p}%)",
"Age (years)" ~ "{median} ({min} to {max})"
)
) %>%
add_p() %>%
bold_labels()
table_IPTW <- survey::svydesign(
id = ~0,
weights = ~weights,
data = data1,
fpc = NULL
) %>%
tbl_svysummary(
by = bs_ab_name,
missing = "no",
statistic = list(
all_continuous() ~ "{median} ({p25}, {p75})",
all_categorical() ~ "{n} ({p}%)"),
include = c("Age (years)", "Female sex","Race","ECOG","LDH (U/L)","Hgb (g/dL)","Baseline ferritin (ng/mL)", "Heavy chain","True EMD","≥1 HRCA","Prior LOTs","Prior ASCT","Prior CAR-T","Triple refractory","Penta refractory","Trial eligible","Prior BCMA-directed therapy", "Received full dose")
) %>%
add_p() %>%
bold_labels()
## The following warnings were returned during `add_p()`:
## ! For variable `Heavy chain` (`bs_ab_name`) and "statistic" and "p.value"
## statistics: Chi-squared approximation may be incorrect
## ! For variable `Race` (`bs_ab_name`) and "statistic" and "p.value" statistics:
## Chi-squared approximation may be incorrect
tbl_merge(list(table_IPTW, table_unweighted), tab_spanner = c("**IPTW**", "**Unweighted**"))
| Characteristic |
IPTW
|
Unweighted
|
||||
|---|---|---|---|---|---|---|
| elra N = 1231 |
tec N = 2481 |
p-value2 | elra N = 1233 |
tec N = 2483 |
p-value4 | |
| Age (years) | 70 (64, 76) | 69 (62, 76) | 0.9 | 71 (39 to 95) | 68 (35 to 88) | 0.060 |
| Female sex | 66 (54%) | 114 (46%) | 0.2 | 68 (55%) | 115 (46%) | 0.11 |
| Race | 0.2 | 0.15 | ||||
| Â Â Â Â White | 94 (79%) | 187 (76%) | 95 (81%) | 187 (76%) | ||
| Â Â Â Â Black | 17 (14%) | 50 (20%) | 16 (14%) | 51 (21%) | ||
| Â Â Â Â AAPI | 3 (2.5%) | 6 (2.3%) | 3 (2.5%) | 5 (2.0%) | ||
| Â Â Â Â American Indian | 2 (1.8%) | 0 (0%) | 2 (1.7%) | 0 (0%) | ||
| Â Â Â Â Other | 2 (2.0%) | 4 (1.8%) | 2 (1.7%) | 4 (1.6%) | ||
| ECOG | >0.9 | 0.5 | ||||
| Â Â Â Â 0 | 23 (18%) | 49 (20%) | 21 (17%) | 50 (20%) | ||
| Â Â Â Â 1 | 62 (50%) | 119 (48%) | 59 (48%) | 121 (49%) | ||
| Â Â Â Â 2 | 28 (23%) | 59 (24%) | 29 (24%) | 60 (24%) | ||
| Â Â Â Â 3+ | 10 (8.2%) | 20 (8.2%) | 14 (11%) | 17 (6.9%) | ||
| LDH (U/L) | 223 (178, 280) | 208 (167, 290) | 0.4 | 222 (175, 277) | 208 (167, 291) | 0.5 |
| Hgb (g/dL) | 9.90 (8.50, 11.40) | 9.80 (8.40, 11.50) | 0.9 | 10.00 (8.70, 11.70) | 9.60 (8.25, 11.20) | 0.027 |
| Baseline ferritin (ng/mL) | 346 (140, 1,250) | 379 (125, 1,219) | >0.9 | 293 (137, 832) | 393 (128, 1,254) | 0.3 |
| Heavy chain | 0.001 | <0.001 | ||||
| Â Â Â Â IgG | 72 (58%) | 109 (44%) | 72 (59%) | 107 (43%) | ||
| Â Â Â Â IgA | 26 (21%) | 43 (17%) | 26 (21%) | 42 (17%) | ||
| Â Â Â Â IgM | 1 (0.7%) | 1 (0.4%) | 1 (0.8%) | 1 (0.4%) | ||
| Â Â Â Â IgD | 2 (1.5%) | 0 (0%) | 2 (1.6%) | 0 (0%) | ||
| Â Â Â Â None | 23 (19%) | 95 (38%) | 22 (18%) | 98 (40%) | ||
| True EMD | 20 (16%) | 39 (16%) | >0.9 | 26 (21%) | 34 (14%) | 0.067 |
| ≥1 HRCA | 81 (70%) | 151 (65%) | 0.3 | 81 (70%) | 152 (65%) | 0.4 |
| Prior LOTs | 6.00 (5.00, 8.00) | 6.00 (4.00, 8.00) | 0.4 | 6.00 (4.00, 8.00) | 6.00 (4.00, 8.00) | >0.9 |
| Prior ASCT | 66 (54%) | 141 (57%) | 0.6 | 61 (50%) | 144 (58%) | 0.12 |
| Prior CAR-T | 53 (43%) | 90 (36%) | 0.2 | 50 (41%) | 94 (38%) | 0.6 |
| Triple refractory | 112 (91%) | 224 (91%) | 0.9 | 112 (91%) | 226 (91%) | >0.9 |
| Penta refractory | 61 (49%) | 123 (49%) | >0.9 | 61 (50%) | 126 (51%) | 0.8 |
| Trial eligible | 24 (19%) | 55 (22%) | 0.5 | 26 (21%) | 53 (21%) | >0.9 |
| Prior BCMA-directed therapy | 67 (54%) | 133 (54%) | >0.9 | 60 (49%) | 138 (56%) | 0.2 |
| Received full dose | 114 (93%) | 241 (97%) | 0.059 | 114 (93%) | 241 (97%) | 0.045 |
| 1 Median (Q1, Q3); n (%) | ||||||
| 2 Design-based KruskalWallis test; Pearson’s X^2: Rao & Scott adjustment | ||||||
| 3 Median (Min to Max); n (%); Median (Q1, Q3) | ||||||
| 4 Wilcoxon rank sum test; Pearson’s Chi-squared test; Fisher’s exact test | ||||||
data1 = IPTW.data %>%
transmute(
BsAb = as.factor(bs_ab_name),
CRS.any = ifelse( max_crs_grade_0_5 != 0, 1, 0),
ICANS.any = ifelse( max_icans_grade != 0,1,0),
CRS.two = ifelse( max_crs_grade_0_5 == 2 | max_crs_grade_0_5 == 3 | max_crs_grade_0_5 == 4, 1, 0),
ICANS.two = ifelse( max_icans_grade ==2 | max_icans_grade == 3 | max_icans_grade == 4,1,0),
CRS.three = ifelse( max_crs_grade_0_5 == 3 | max_crs_grade_0_5 == 4, 1, 0),
ICANS.three = ifelse( max_icans_grade == 3 | max_icans_grade == 4,1,0),
weights
)
weighted_table <- survey::svydesign(
id = ~0,
weights = ~weights,
data = data1,
fpc = NULL
)
combined_data <- rbind(
IPTW.data %>%
count(bs_ab_name, max_crs_grade_0_5, wt=weights) %>%
group_by(bs_ab_name) %>%
mutate(
bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
total = sum(n[max_crs_grade_0_5 != 0]),
pct = prop.table(n) * 100,
total_pct = sum(pct[max_crs_grade_0_5 != 0]),
Grade = factor(max_crs_grade_0_5, levels = c(0, 4, 3, 2, 1)),
Category = "CRS"
),
IPTW.data %>%
count(bs_ab_name, max_icans_grade, wt=weights) %>%
group_by(bs_ab_name) %>%
mutate(
bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
total = sum(n[max_icans_grade != 0]),
pct = prop.table(n) * 100,
total_pct = sum(pct[max_icans_grade != 0]),
Grade = factor(max_icans_grade, levels = c(0, 4, 3, 2, 1)),
Category = "ICANS"
)
) %>%
mutate(
Category = factor(Category, levels = c("CRS","ICANS"))
)
# Any grade toxicity
CRS.pval <- as.numeric(svychisq(~CRS.any+BsAb, weighted_table)$p.value)
ICANS.pval <- as.numeric(svychisq(~ICANS.any+BsAb, weighted_table)$p.value)
f_labels <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.pval, ICANS.pval) ) %>%
mutate(
Category = factor(Category, levels = c("CRS","ICANS"))
)
# G2+ toxicity
CRS.two.pval <- as.numeric(svychisq(~CRS.two+BsAb, weighted_table)$p.value)
ICANS.two.pval <- as.numeric(svychisq(~ICANS.two+BsAb, weighted_table)$p.value)
f_labels_two <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.two.pval, ICANS.two.pval) ) %>%
mutate(
Category = factor(Category, levels = c("CRS","ICANS"))
)
# G3+ toxicity
CRS.three.pval <- as.numeric(svychisq(~CRS.three+BsAb, weighted_table)$p.value)
ICANS.three.pval <- as.numeric(svychisq(~ICANS.three+BsAb, weighted_table)$p.value)
f_labels_three <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.three.pval, ICANS.three.pval) ) %>%
mutate(
Category = factor(Category, levels = c("CRS","ICANS"))
)
toxplot_IPTW <- ggplot(combined_data) +
aes(x = bs_ab_name, y = pct, fill = Grade) +
geom_bar(stat = "identity", position = "stack", data = . %>% filter(Grade != 0)) +
geom_text(aes(label = paste0(signif(pct,2), "% (", round(n,0) ,")")),
data = . %>% filter(Grade != 0 & Grade != 4),
position = position_stack(vjust = 0.5)) +
geom_text( aes(label = paste0(signif(total_pct,2), "% (", signif(total,2),")"), x = bs_ab_name, y = total_pct + 6, vjust = 0),
data = . %>% filter(Grade == 1),
color = "#104a8e",
fontface = "bold") +
facet_wrap(~Category, scales = "free_y") +
ylab("Proportion of patients (% [n])") +
xlab(NULL) +
theme_classic()+
scale_fill_brewer(palette = "Pastel1") +
#guides(fill = guide_legend(reverse = TRUE)) +
scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 110)) +
theme(
plot.title = element_text(size = 14),
axis.title = element_text(size = 14),
axis.text = element_text(size = 14),
legend.title = element_text(size = 14),
legend.text = element_text(size = 14),
strip.text = element_text(size = 14),
axis.text.x = element_text(angle = 45, hjust = 1),
#legend.position = "bottom"
legend.position = "right"
) +
coord_cartesian(ylim = c(0, 115)) +
geom_text( aes(label = paste0("Grade 1+, p=", signif(label,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels ) +
geom_text( aes(label = paste0("Grade 2+, p=", signif(label,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_two ) +
geom_text( aes(label = paste0("Grade 3+, p=", signif(label,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_three ) +
ggtitle("A) IPTW")
toxplot_IPTW
combined_data <- rbind(
combined_db %>%
count(bs_ab_name, max_crs_grade_0_5) %>%
group_by(bs_ab_name) %>%
mutate(
total = sum(n[max_crs_grade_0_5 != 0]),
pct = prop.table(n) * 100,
total_pct = sum(pct[max_crs_grade_0_5 != 0]),
Grade = factor(max_crs_grade_0_5, levels = c(0, 4, 3, 2, 1)),
bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
Category = "CRS"
),
combined_db %>%
count(bs_ab_name, max_icans_grade) %>%
group_by(bs_ab_name) %>%
mutate(
total = sum(n[max_icans_grade != 0]),
pct = prop.table(n) * 100,
total_pct = sum(pct[max_icans_grade != 0]),
Grade = factor(max_icans_grade, levels = c(0, 4, 3, 2, 1)),
bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
Category = "ICANS"
)
) %>%
mutate(
Category = factor(Category, levels = c("CRS","ICANS"))
)
## Any-grade toxicity
CRS.d.test <- data.frame(
results = ifelse(combined_db$max_crs_grade_0_5 !=0,1,0) ,
test = combined_db$bs_ab_name
)
CRS.pval <- chisq.test(CRS.d.test$results,CRS.d.test$test, correct = FALSE)$p.value
ICANS.d.test <- data.frame(
results = ifelse(combined_db$max_icans_grade !=0,1,0) ,
test = combined_db$bs_ab_name
)
ICANS.pval <- chisq.test(ICANS.d.test$results,ICANS.d.test$test, correct = FALSE)$p.value
f_labels <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.pval, ICANS.pval) ) %>%
mutate(
Category = factor(Category, levels = c("CRS","ICANS") )
)
## G2+ toxicity
CRS.two.d.test <- data.frame(
results = ifelse(combined_db$max_crs_grade_0_5 == 2 | combined_db$max_crs_grade_0_5 == 3 | combined_db$max_crs_grade_0_5 == 4,1,0),
test = combined_db$bs_ab_name
)
CRS.two.pval <- chisq.test(CRS.two.d.test$results,CRS.two.d.test$test, correct = FALSE)$p.value
ICANS.two.d.test <- data.frame(
results = ifelse(combined_db$max_icans_grade == 2 | combined_db$max_icans_grade == 3 | combined_db$max_icans_grade == 4,1,0),
test = combined_db$bs_ab_name
)
ICANS.two.pval <- chisq.test(ICANS.two.d.test$results,ICANS.two.d.test$test, correct = FALSE)$p.value
f_labels_two <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.two.pval, ICANS.two.pval) ) %>%
mutate(
Category = factor(Category, levels = c("CRS","ICANS") )
)
## G3+ toxicity
CRS.three.d.test <- data.frame(
results = ifelse( combined_db$max_crs_grade_0_5 == 3 | combined_db$max_crs_grade_0_5 == 4,1,0),
test = combined_db$bs_ab_name
)
CRS.three.pval <- fisher.test(CRS.three.d.test$results,CRS.three.d.test$test)$p.value
ICANS.three.d.test <- data.frame(
results = ifelse(combined_db$max_icans_grade == 3 | combined_db$max_icans_grade == 4,1,0),
test = combined_db$bs_ab_name
)
ICANS.three.pval <- chisq.test(ICANS.three.d.test$results,ICANS.three.d.test$test, correct = FALSE)$p.value
f_labels_three <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.three.pval, ICANS.three.pval) ) %>%
mutate(
Category = factor(Category, levels = c("CRS","ICANS") )
)
toxplot_unweighted <- ggplot(combined_data) +
aes(x = bs_ab_name, y = pct, fill = Grade) +
geom_bar(stat = "identity", position = "stack", data = . %>% filter(Grade != 0)) +
geom_text(aes(label = paste0(signif(pct,2), "% (",n,")")),
data = . %>% filter(Grade != 0 & Grade != 4),
position = position_stack(vjust = 0.5)) +
geom_text( aes(label = paste0(signif(total_pct,2), "% (",total,")"), x = bs_ab_name, y = total_pct + 6, vjust = 0),
data = . %>% filter(Grade == 1),
color = "#104a8e",
fontface = "bold") +
facet_wrap(~Category, scales = "free_y") +
ylab("Proportion of patients (% [n])") +
xlab(NULL) +
theme_classic()+
scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 100)) +
theme(
plot.title = element_text(size = 14),
axis.title = element_text(size = 14),
axis.text = element_text(size = 14),
legend.title = element_text(size = 14),
legend.text = element_text(size = 14),
strip.text = element_text(size = 14),
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none"
) +
coord_cartesian(ylim = c(0, 114)) +
scale_fill_brewer(palette = "Pastel1", direction = 1) +
geom_text( aes(label = paste0("Grade 1+, p=", signif(label,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels ) +
geom_text( aes(label = paste0("Grade 2+, p=", signif(label,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_two ) +
geom_text( aes(label = paste0("Grade 3+, p=", signif(label,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_three ) +
ggtitle("B) Unweighted")
toxplot_unweighted
FIGURE1 <- toxplot_IPTW / toxplot_unweighted
ggsave("svg/FIGURE1.svg", FIGURE1, device = "svg",
width = 7,
height = 11,
units = c("in"))
ggsave("tiff/FIGURE1.tiff", FIGURE1, device = "tiff",
width = 7,
height = 11,
dpi = 600,
units = c("in"),
compression = "lzw")
ggsave("jpeg/FIGURE1.jpeg", FIGURE1, device = "jpeg",
width = 7,
height = 11,
dpi = 600,
units = "in")
knitr::include_graphics("jpeg/FIGURE1.jpeg")
data1 = IPTW.data %>%
filter(!is.na(best_response)) %>%
transmute(
BsAb = as.factor(bs_ab_name),
best_response = factor(best_response, levels = c("sCR", "CR","VGPR","PR","MR","SD","PD")),
ORR = ifelse(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR",1,0),
VGPR = ifelse(best_response == "sCR" | best_response == "CR" | best_response == "VGPR",1,0),
CR = ifelse(best_response == "sCR" | best_response == "CR",1,0),
weights
)
weighted_table <- survey::svydesign(
id = ~0,
weights = ~weights,
data = data1,
fpc = NULL
)
data1 <- IPTW.data %>%
filter(!is.na(best_response)) %>%
count(bs_ab_name, best_response, wt=weights) %>%
group_by(bs_ab_name) %>%
mutate(
total = sum(n[best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"]),
pct = prop.table(n) * 100,
total_pct = sum(pct[best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"]),
best_response = factor(best_response, levels = c("sCR", "CR","VGPR","PR","MR", "SD","PD")),
bs_ab_name = as.factor(bs_ab_name),
Category = "Response rate"
)
Response.pval <- as.numeric(svychisq(~ORR+BsAb, weighted_table)$p.value)
CR.pval <- as.numeric(svychisq(~CR+BsAb, weighted_table)$p.value)
VGPR.pval <- as.numeric(svychisq(~VGPR+BsAb, weighted_table)$p.value)
f_labels <- data.frame(Category = c("Response rate"), label = c(Response.pval) ) %>%
mutate(
Category = factor(Category, levels = c("Response rate"))
)
f_labels_CR <- data.frame(Category = c("Response rate"), label = c(CR.pval) ) %>%
mutate(
Category = factor(Category, levels = c("Response rate"))
)
f_labels_VGPR <- data.frame(Category = c("Response rate"), label = c(VGPR.pval) ) %>%
mutate(
Category = factor(Category, levels = c("Response rate"))
)
responseplot_IPTW <- ggplot(data1) +
aes(x = bs_ab_name, y = pct, fill = best_response) +
geom_bar(stat = "identity", position = "stack", data = . %>% filter(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR")) +
geom_text(aes(label = paste0(signif(pct,2), "% (", round(n,0) ,")")),
data = . %>% filter(best_response == "sCR" | best_response == "CR" | best_response == "VGPR"| best_response == "PR"),
position = position_stack(vjust = 0.5)) +
geom_text( aes(label = paste0(signif(total_pct,2), "% (", signif(total,2),")"), x = bs_ab_name, y = total_pct + 7, vjust = 0),
data = . %>% filter(best_response == "PR"),
color = "#104a8e",
fontface = "bold") +
facet_wrap(~Category, scales = "free_y") +
ylab("Proportion of patients (% [n])") +
xlab(NULL) +
theme_classic()+
scale_fill_brewer(
palette = "Pastel1",
direction = -1,
breaks = rev(levels(data1$best_response))
) +
guides(fill = guide_legend(reverse = TRUE)) +
scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 110)) +
theme(
plot.title = element_text(size = 14),
axis.title = element_text(size = 14),
axis.text = element_text(size = 14),
legend.title = element_text(size = 14),
legend.text = element_text(size = 14),
strip.text = element_text(size = 14),
axis.text.x = element_text(angle = 45, hjust = 1),
#legend.position = "bottom"
legend.position = "right"
) +
coord_cartesian(ylim = c(0, 114)) +
geom_text( aes(label = paste0("ORR, p=", signif(Response.pval,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels ) +
geom_text( aes(label = paste0("≥VGPR, p=", signif(VGPR.pval,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_VGPR ) +
geom_text( aes(label = paste0("≥CR, p=", signif(CR.pval,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_CR ) +
labs(fill = "Best response")+
ggtitle("A) IPTW")
responseplot_IPTW
data1 <- combined_db %>%
filter(!is.na(best_response)) %>%
transmute(
best_response = factor(
best_response,
levels = c("sCR", "CR","VGPR", "PR","MR","SD","PD")),
ORR = ifelse(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" |best_response == "PR",1,0),
bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
Category = "Response rate"
)
combined_data <- combined_db %>%
filter(!is.na(best_response)) %>%
count(bs_ab_name, best_response) %>%
group_by(bs_ab_name) %>%
transmute(
bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
best_response = factor(
best_response,
levels = c("sCR", "CR","VGPR", "PR","MR","SD","PD")),
total = sum(n[best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"]),
pct = prop.table(n) * 100,
n,
total_pct = sum(pct[best_response == "sCR" | best_response == "CR" |best_response == "VGPR" | best_response == "PR"]),
Category = "Response rate"
)
## ORR
Response.d.test <- data.frame(
results = ifelse(data1$ORR,1,0),
test = data1$bs_ab_name
)
Response.pval <- chisq.test(Response.d.test$results,Response.d.test$test, correct = FALSE)$p.value
## VGPR
VGPR.d.test <- data.frame(
results = ifelse(data1$best_response == "CR" | data1$best_response == "VGPR",1,0),
test = data1$bs_ab_name
)
VGPR.pval <- chisq.test(VGPR.d.test$results,VGPR.d.test$test, correct = FALSE)$p.value
## CR
CR.d.test <- data.frame(
results = ifelse(data1$best_response == "CR",1,0),
test = data1$bs_ab_name
)
CR.pval <- chisq.test(CR.d.test$results,CR.d.test$test, correct = FALSE)$p.value
f_labels <- data.frame(Category = c("Response rate"), label = c(Response.pval) ) %>%
mutate(
Category = factor(Category, levels = c("Response rate"))
)
f_labels_CR <- data.frame(Category = c("Response rate"), label = c(CR.pval) ) %>%
mutate(
Category = factor(Category, levels = c("Response rate"))
)
f_labels_VGPR <- data.frame(Category = c("Response rate"), label = c(VGPR.pval) ) %>%
mutate(
Category = factor(Category, levels = c("Response rate"))
)
responseplot_unweighted <- ggplot(combined_data) +
aes(x = bs_ab_name, y = pct, fill = best_response) +
geom_bar(stat = "identity", position = "stack", data = . %>% filter(best_response == "sCR" | best_response == "CR" |best_response == "VGPR" | best_response == "PR")) +
geom_text(aes(label = paste0(signif(pct,2), "% (",n,")")),
data = . %>% filter(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"),
position = position_stack(vjust = 0.5)) +
geom_text( aes(label = paste0(signif(total_pct,2), "% (",total,")"), x = bs_ab_name, y = total_pct + 7, vjust = 0),
data = . %>% filter(best_response == "PR"),
color = "#104a8e",
fontface = "bold") +
facet_wrap(~Category, scales = "free_y") +
ylab("Proportion of patients (% [n])") +
xlab(NULL) +
theme_classic()+
scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 110)) +
theme(
plot.title = element_text(size = 14),
axis.title = element_text(size = 14),
axis.text = element_text(size = 14),
legend.title = element_text(size = 14),
legend.text = element_text(size = 14),
strip.text = element_text(size = 14),
axis.text.x = element_text(angle = 45, hjust = 1),
#legend.position = "bottom"
legend.position = "none"
) +
coord_cartesian(ylim = c(0, 114)) +
scale_fill_brewer(palette = "Pastel1", direction = -1) +
geom_text( aes(label = paste0("ORR, p=", signif(Response.pval,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels) +
geom_text( aes(label = paste0("≥VGPR, p=", signif(VGPR.pval,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_VGPR ) +
geom_text( aes(label = paste0("≥CR, p=", signif(CR.pval,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_CR ) +
labs(fill = "Best response")+
ggtitle("B) Unweighted")
responseplot_unweighted
FIGURE2 <- responseplot_IPTW / responseplot_unweighted
ggsave("svg/FIGURE2.svg", FIGURE2, device = "svg",
width = 7,
height = 11,
units = c("in"))
ggsave("tiff/FIGURE2.tiff", FIGURE2, device = "tiff",
width = 7,
height = 11,
dpi = 600,
units = c("in"),
compression = "lzw")
ggsave("jpeg/FIGURE2.jpeg", FIGURE2, device = "jpeg",
width = 7,
height = 11,
dpi = 600,
units = "in")
knitr::include_graphics("jpeg/FIGURE2.jpeg")
data1 <- IPTW.data %>%
transmute(
bs_ab_name,
Death = ifelse(death_yes_no,1,0),
Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
Reverse_death = ifelse(death_yes_no == 1, 0,1),
weights
)
reverse_km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 1, weights = weights, data1)
reverse_km_OS
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Reverse_death) ~
## 1, data = data1, weights = weights)
##
## 15 observations deleted due to missingness
## records n events median 0.95LCL 0.95UCL
## [1,] 356 356 215 13.4 12.5 14.9
reverse_km_OS_elra <- survfit(Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 1, weights = weights, data1 %>% filter(bs_ab_name == "elra"))
reverse_km_OS_elra
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Reverse_death) ~
## 1, data = data1 %>% filter(bs_ab_name == "elra"), weights = weights)
##
## 6 observations deleted due to missingness
## records n events median 0.95LCL 0.95UCL
## [1,] 117 117 74.5 7.47 6.03 9
reverse_km_OS_tec <- survfit(Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 1, weights = weights, data1 %>% filter(bs_ab_name == "tec"))
reverse_km_OS_tec
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Reverse_death) ~
## 1, data = data1 %>% filter(bs_ab_name == "tec"), weights = weights)
##
## 9 observations deleted due to missingness
## records n events median 0.95LCL 0.95UCL
## [1,] 239 239 140 17.1 14.9 20.1
data1 <- IPTW.data %>%
transmute(
site,
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
bs_ab_name,
Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
Death = ifelse(death_yes_no,1,0),
date_of_bs_ab_step_up_d1,
date_of_pd,
date_of_last_contact,
date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death - date_of_bs_ab_step_up_d1),
trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
weights
)
data_dor <- IPTW.data %>%
filter(best_response %in% c("sCR","CR","VGPR","PR")) %>%
transmute(
weights,
bs_ab_name,
day_30_response,
day_90_response,
best_response,
Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
Death = ifelse(death_yes_no,1,0),
date_of_pd,
date_of_last_contact,
date_of_bs_ab_step_up_d1,
date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
day_30_response,
day_90_response,
best_response,
time_to_best_response = case_when(
best_response == day_30_response ~ 30,
best_response == day_90_response ~ 90,
.default = NA
),
trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0)
) %>%
filter(!is.na(time_to_best_response)) %>%
mutate(
Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death - date_of_bs_ab_step_up_d1)- time_to_best_response,
) %>%
filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)
km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, weights = weights, data = data1)
km_OS
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name,
## data = data1, weights = weights)
##
## 15 observations deleted due to missingness
## records n events median 0.95LCL 0.95UCL
## bs_ab_name=elra 117 117 42.8 NA 8.6 NA
## bs_ab_name=tec 239 239 98.6 21.5 16.4 NA
med_time <- surv_median(km_OS)
cox_model <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)
HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)
## 1-year timepoint
timepoint1 <- summary(km_OS,times=c(3))
timepoint2 <- summary(km_OS,times=c(6))
timepoint3 <- summary(km_OS,times=c(9))
timepoint4 <- summary(km_OS,times=c(12))
timepoint5 <- summary(km_OS,times=c(15))
OS <- ggsurvplot(km_OS,
pval=FALSE,
conf.int = FALSE,
risk.table = TRUE, # Add risk table
fontsize = 4,
risk.table.col = "strata", # Change risk table color by groups
tables.height = 0.3,
tables.theme = theme_cleantable() +
theme(
axis.text.y = element_text(size = 9),
plot.title = element_text(size = 10)
),
xlab="Time (months)", ylab = "Overall survival",
xlim = c(0,18), break.x.by = c(3),
ylim = c(0,1), break.y.by = c(0.25),
legend = "none",
surv.scale="percent",
font.main = c(10, "plain", "black"),
font.x = c(10, "plain", "black"),
font.y = c(10, "plain", "black"),
font.caption = c(10, "plain", "black"),
font.tickslab = c(10, "plain", "black"),
legend.labs = c("Elra", "Tec"),
title = "A)",
size = 0.5
)
OS$plot <- OS$plot +
geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
annotate("text", x = 3, y = timepoint1$surv[1] - 0.1, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = -0.7, y = 0,
label = paste0("Elra vs tec\n",
" Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
" HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
color = "black", hjust = 0, vjust = 0, size = 3)
OS
km_PFS <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data1)
km_PFS
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~
## bs_ab_name, data = data1, weights = weights)
##
## 12 observations deleted due to missingness
## records n events median 0.95LCL 0.95UCL
## bs_ab_name=elra 118 118 68.1 3.77 2.57 7.17
## bs_ab_name=tec 241 241 151.1 8.50 7.10 11.30
med_time <- surv_median(km_PFS)
cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)
HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)
## Timepoints
timepoint1 <- summary(km_PFS,times=c(3))
timepoint2 <- summary(km_PFS,times=c(6))
timepoint3 <- summary(km_PFS,times=c(9))
timepoint4 <- summary(km_PFS,times=c(12))
timepoint5 <- summary(km_PFS,times=c(15))
PFS <- ggsurvplot(km_PFS,
pval=FALSE,
conf.int = FALSE,
risk.table = TRUE, # Add risk table
fontsize = 4,
risk.table.col = "strata", # Change risk table color by groups
tables.height = 0.3,
tables.theme = theme_cleantable() +
theme(
axis.text.y = element_text(size = 9),
plot.title = element_text(size = 10)
),
xlab="Time (months)", ylab = "Progression-free survival",
xlim = c(0,18), break.x.by = c(3),
ylim = c(0,1), break.y.by = c(0.25),
legend = "none",
surv.scale="percent",
font.main = c(10, "plain", "black"),
font.x = c(10, "plain", "black"),
font.y = c(10, "plain", "black"),
font.caption = c(10, "plain", "black"),
font.tickslab = c(10, "plain", "black"),
legend.labs = c("Elra", "Tec"),
title = "B)",
size = 0.5
)
PFS$plot <- PFS$plot +
geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
annotate("text", x = 3, y = timepoint1$surv[1] - 0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = -0.7, y = 0,
label = paste0("Elra vs tec\n",
" Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
" HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
color = "black", hjust = 0, vjust = 0, size = 3)
PFS
km_DOR <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data_dor)
km_DOR
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~
## bs_ab_name, data = data_dor, weights = weights)
##
## 3 observations deleted due to missingness
## records n events median 0.95LCL 0.95UCL
## bs_ab_name=elra 62 60.2 20.6 11.7 6.13 NA
## bs_ab_name=tec 118 120.7 56.1 13.7 9.27 22.3
med_time <- surv_median(km_DOR)
cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data_dor)
cox_summary <- summary(cox_model)
HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)
## Specific timepoints
timepoint1 <- summary(km_DOR,times=c(3))
timepoint2 <- summary(km_DOR,times=c(6))
timepoint3 <- summary(km_DOR,times=c(9))
timepoint4 <- summary(km_DOR,times=c(12))
timepoint5 <- summary(km_DOR,times=c(15))
DOR <- ggsurvplot(km_DOR,
pval=FALSE,
conf.int = FALSE,
risk.table = TRUE, # Add risk table
fontsize = 4,
risk.table.col = "strata", # Change risk table color by groups
tables.height = 0.3,
#risk.table.fontsize = 4,
tables.theme = theme_cleantable() +
theme(
axis.text.y = element_text(size = 9),
plot.title = element_text(size = 10)
),
xlab="Time (months)", ylab = "Duration of response",
xlim = c(0,18), break.x.by = c(3),
ylim = c(0,1), break.y.by = c(0.25),
legend = "none",
surv.scale="percent",
font.main = c(10, "plain", "black"),
font.x = c(10, "plain", "black"),
font.y = c(10, "plain", "black"),
font.caption = c(10, "plain", "black"),
font.tickslab = c(10, "plain", "black"),
legend.labs = c("Elra", "Tec"),
title = "C)",
size = 0.5
)
DOR$plot <- DOR$plot +
geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
annotate("text", x = 3, y = timepoint1$surv[1] -0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = -0.7, y = 0,
label = paste0("Elra vs tec\n",
" Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
" HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
color = "black", hjust = 0, vjust = 0, size = 3)
DOR
FIGURE3 <- wrap_plots(
(OS$plot / OS$table) + plot_layout(heights = c(4, 0.9)),
(PFS$plot / PFS$table) + plot_layout(heights = c(4, 0.9)),
(DOR$plot / DOR$table) + plot_layout(heights = c(4, 0.9)),
ncol = 1 # Force vertical stacking
)
ggsave("svg/FIGURE3.svg", FIGURE3, device = "svg",
width = 6.5,
height = 8.5,
units = c("in"))
ggsave("jpeg/FIGURE3.jpeg", FIGURE3, device = "jpeg",
width = 6.5,
height = 8.5,
dpi = 600,
units = c("in"))
ggsave("tiff/FIGURE3.tiff", FIGURE3, device = "tiff",
width = 6.5,
height = 8.5,
dpi = 600,
units = c("in"),
compression = "lzw")
knitr::include_graphics("jpeg/FIGURE3.jpeg")
data1 <- IPTW.data %>%
transmute(
site,
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
elra_ind = ifelse(bs_ab_name == "elra", 1, 0), # 1 = elra, 0 = tec
Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
Death = ifelse(death_yes_no,1,0),
date_of_bs_ab_step_up_d1,
date_of_pd,
date_of_last_contact,
date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death - date_of_bs_ab_step_up_d1),
trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
weights
)
## ---- Cox models ----
cox_OS <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ elra_ind, weights = weights, data = data1)
cox_PFS <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ elra_ind, weights = weights, data = data1)
## ---- PH diagnostics ----
zph_OS <- cox.zph(cox_OS)
zph_PFS <- cox.zph(cox_PFS)
# Global PH p-values (same as elra_ind right now)
p_OS <- signif(zph_OS$table["GLOBAL", "p"], 3)
p_PFS <- signif(zph_PFS$table["GLOBAL", "p"], 3)
par(mfrow = c(1, 2))
plot(
zph_OS,
main = paste0(
"Schoenfeld residuals – OS\n",
"Global PH test: p = ", p_OS
)
)
plot(
zph_PFS,
main = paste0(
"Schoenfeld residuals – PFS\n",
"Global PH test: p = ", p_PFS
)
)
par(mfrow = c(1, 1))